# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import warnings
from hysop.topology.topology import (
Topology,
TopologyState,
TopologyView,
TopologyWarning,
)
from hysop.constants import np, math, Backend, MemoryOrdering
from hysop.constants import HYSOP_ORDER, BoundaryCondition, HYSOP_INTEGER
from hysop.constants import HYSOP_REAL, HYSOP_MPI_REAL, TranspositionState
from hysop.tools.transposition_states import TranspositionState
from hysop.domain.box import Box, BoxView
from hysop.core.mpi import MPI
from hysop.tools.htypes import check_instance, to_tuple, first_not_None
from hysop.tools.parameters import CartesianDiscretization, MPIParams
from hysop.tools.misc import Utils, prod
from hysop.tools.decorators import debug, deprecated
from hysop.tools.numpywrappers import npw
from hysop.tools.string_utils import prepend
from hysop.mesh.cartesian_mesh import CartesianMesh, CartesianMeshView
from hysop.testsenv import __HAS_OPENCL_BACKEND__
[docs]
class CartesianTopologyState(TopologyState):
"""
CartesianTopology topology state.
This is a helper class to qualify CartesianDiscreteField states.
A CartesianTopologyState contains informations about
the way the application should perceive the contained data
(Arrays) in CartesianDiscreteFields.
Those informations include for the current physical
transposition state of the topology and the local meshes
and the memory_order.
Currently the state is shared accross all components of the
CartesianDiscreteField. This is global state for all processes
contained in the linked Cartesian topology.
"""
__slots__ = ("_is_read_only", "_dim", "_axes", "_memory_order", "_task_id")
@debug
def __new__(
cls, dim, task_id, axes=None, memory_order=None, is_read_only=False, **kwds
):
return super().__new__(cls, is_read_only=is_read_only, **kwds)
@debug
def __init__(
self, dim, task_id, axes=None, memory_order=None, is_read_only=False, **kwds
):
"""
Initialize a CartesianState to given parameters.
Parameters
----------
dim: int
Dimensiong of the underlying domain.
axes: tuple of ints, optional
Permutation of the CartesianTopology and the meshes as tuple of ints,
in numpy notations.
memory_order: :class:`hysop.constants.MemoryOrdering`, optional
Either C or Fortran memory order.
is read_only: bool, optional
Whether this topology is read only or not.
kwds: dict
Base class keyword arguments.
"""
super().__init__(is_read_only=is_read_only, **kwds)
self._dim = int(dim)
self._task_id = int(task_id)
self._set_axes(axes)
self._set_memory_order(memory_order)
def _get_axes(self):
"""Return current permutation as a tuple of int."""
return self._axes
def _set_axes(self, axes):
"""Set the current permutation as a tuple of int."""
axes = axes or TranspositionState[self._dim].default_axes()
axes = to_tuple(axes, cast=int)
check_instance(axes, tuple, values=int)
assert set(axes) == set(range(self._dim))
self._axes = axes
def _get_memory_order(self):
"""Return current memory_order as a tuple (one memory_order per axis)."""
return self._memory_order
def _set_memory_order(self, memory_order):
"""Set the current memory_order as a tuple of hysop.constants.memory_order."""
memory_order = first_not_None(memory_order, MemoryOrdering.C_CONTIGUOUS)
if memory_order == "c":
memory_order = MemoryOrdering.C_CONTIGUOUS
elif memory_order == "f":
memory_order = MemoryOrdering.F_CONTIGUOUS
assert memory_order in (
MemoryOrdering.C_CONTIGUOUS,
MemoryOrdering.F_CONTIGUOUS,
), memory_order
self._memory_order = memory_order
def _get_dim(self):
"""Return the dimension of the underlying topology domain."""
return self._dim
def _get_task_id(self):
"""Return the task identifier of the underlying topology domain."""
return self._task_id
def _get_tstate(self):
"""Return the TranspositionState corresponding to current permutation axes."""
return TranspositionState.axes_to_tstate(self._axes)
def __transposed(self, vec, axes):
"""
Compute permutation of input vector of size len(axes) according
to axes permutation.
"""
axes = to_tuple(axes)
if self.memory_order is MemoryOrdering.F_CONTIGUOUS:
axes = axes[::-1]
assert len(axes) > 0
assert set(axes) == set(range(len(axes))), axes
if isinstance(vec, np.ndarray):
if (vec.ndim > 1) and (vec.ndim == len(axes)):
assert vec.size != len(axes), "ambiguous transposition"
return npw.transpose(vec, axes=axes)
else:
assert vec.size == len(axes), f"{vec.size} != {len(axes)}"
res = vec.copy()
res[...] = tuple(vec[i] for i in axes)
return res
else:
assert len(vec) == len(axes)
return type(vec)(vec[i] for i in axes)
[docs]
def transposed(self, vec):
"""Compute permutation of input vector according to current transposition state."""
return self.__transposed(vec, self._axes)
[docs]
def copy(self, axes=None, memory_order=None, is_read_only=None):
"""Return a copy of self, some properties may be alterted in kwds."""
memory_order = first_not_None(memory_order, self.memory_order)
is_read_only = first_not_None(is_read_only, self.is_read_only)
axes = first_not_None(axes, self.axes)
return CartesianTopologyState(
dim=self.dim,
task_id=self.task_id,
axes=axes,
memory_order=memory_order,
is_read_only=is_read_only,
)
[docs]
def short_description(self):
"""Return a short description of this CartesianTopologyState."""
s = "{}[order={}, axes=({}), ro={}, task={}]"
return s.format(
self.full_tag,
self.memory_order,
",".join(str(a) for a in self.axes),
"1" if self.is_read_only else "0",
self.task_id,
)
[docs]
def long_description(self):
"""Return a long description of this CartesianTopologyState."""
s = """{}
*dim: {}
*task_id: {}
*order: {}
*axes: ({})
*tstate: {}
*read only: {}
"""
return s.format(
self.full_tag,
self.dim,
self.task_id,
self.memory_order,
",".join([str(a) for a in self.axes]),
self.tstate,
self.is_read_only,
)
[docs]
def match(self, other, invert=False):
"""Check if this topology state does match the other one."""
if not isinstance(other, CartesianTopologyState):
return NotImplemented
match = super().match(other, invert=False)
match &= self._dim == other._dim
match &= self._task_id == other._task_id
match &= self._axes == other._axes
match &= self._memory_order == other._memory_order
return not match if invert else match
def __hash__(self):
h = super().__hash__()
return (
h
^ hash(self._dim)
^ hash(self._task_id)
^ hash(self._axes)
^ hash(self._memory_order)
)
dim = property(_get_dim)
task_id = property(_get_task_id)
tstate = property(_get_tstate)
axes = property(_get_axes, _set_axes)
memory_order = property(_get_memory_order, _set_memory_order)
[docs]
class CartesianTopologyView(TopologyView):
"""
CartesianTopology topology view.
Basically a CartesianTopology topology with field transposition state taken into account.
"""
__slots__ = ("_mesh_view", "_domain_view", "_topology", "_topology_state")
@debug
def __init__(self, topology_state, topology=None, **kwds):
super().__init__(topology_state=topology_state, topology=topology, **kwds)
[docs]
@debug
def __new__(cls, topology_state, topology=None, **kwds):
"""
Create and initialize a cartesian topology view.
Parameters
----------
topology_state: :class:`~hysop.topology.cartesian_topology.CartesianTopologyState`
State that charaterizes the given view.
topology: :class:`~hysop.topology.topology.CartesianTopology`
Original cartesian topology on which the view is.
kwds: dict
Base class keyword arguments.
"""
check_instance(topology_state, CartesianTopologyState)
check_instance(topology, CartesianTopology, allow_none=True)
obj = super().__new__(
cls, topology_state=topology_state, topology=topology, **kwds
)
check_instance(obj._topology, CartesianTopology)
return obj
def _get_proc_axes(self):
"""Returns state transposition axes."""
return self._topology_state.axes
def _get_proc_memory_order(self):
"""Returns state memory_order."""
return self._topology_state.memory_order
def _proc_transposed(self, vec):
"""Returns transposed vec according to current transposition state."""
return self._topology_state.transposed(vec)
def _distributed_components(self, vec):
"""
Extract distributed components of vector and returns
reduced vector of same type.
"""
rvec = np.asarray(vec)[self._get_is_distributed()]
if isinstance(vec, np.ndarray):
return rvec.copy()
else:
return type(vec)(rvec.tolist())
# ATTRIBUTE GETTERS
def _get_global_resolution(self):
"""Returns global resolution of the discretization (logical grid size)."""
return self._proc_transposed(self._topology._discretization.global_resolution)
def _get_grid_resolution(self):
"""Returns grid resolution of the discretization (effective grid size)."""
return self._proc_transposed(self._topology._discretization.grid_resolution)
def _get_ghosts(self):
"""Returns ghosts of the discretization."""
return self._proc_transposed(self._topology._discretization.ghosts)
def _get_comm(self):
"""Return the communicator of this topology."""
return self._get_cart_comm()
def _get_cart_comm(self):
"""MPI cartesian topology intracommunicator."""
return self._topology._cart_comm
def _get_cart_dim(self):
"""MPI cartesian topology dimension."""
return self._topology._cart_dim
def _get_cart_size(self):
"""MPI cartesian topology size."""
return self._topology._cart_size
def _get_cart_shape(self):
"""MPI cartesian topology shape."""
return self._distributed_components(self._get_proc_shape())
def _get_cart_periods(self):
"""MPI cartesian topology shape."""
return self._distributed_components(self._get_is_periodic())
def _get_cart_coords(self):
"""Current process MPI cartesian topology coordinates."""
return self._distributed_components(self._get_proc_coords())
def _get_cart_rank(self):
"""Current process MPI cartesian topology rank."""
return self._topology._cart_rank
def _get_cart_ranks(self):
"""Return all MPI cartesian topology ranks as np.ndarray."""
return self._get_proc_ranks().reshape(self._get_cart_shape())
def _get_cart_ranks_mapping(self):
"""Return all parent communicator topology ranks as np.ndarray."""
return self._get_proc_ranks_mapping().reshape(self._get_cart_shape())
def _get_cart_neighbour_ranks(self):
"""Return the ranks of the neighbours nodes as obtained by MPI_Cart_shift.
self.neighbours[0,i] (resp. [1,i]) is the previous (resp. next)
neighbour in axe i."""
return self._get_proc_neighbour_ranks()[:, self._get_is_distributed()]
def _get_proc_shape(self):
"""MPI cartesian topology extended shape (ie. undistributed axes included)."""
return self._proc_transposed(self._topology._proc_shape)
def _get_proc_coords(self):
"""
Current process cartesian topology extended coordinates (ie. undistributed
axes included).
"""
return self._proc_transposed(self._topology._proc_coords)
def _get_proc_ranks(self):
"""
Return all MPI cartesian topology ranks as np.ndarray
undistributed axes included.
"""
return np.transpose(self._topology._proc_ranks, axes=self._get_proc_axes())
def _get_proc_ranks_mapping(self):
"""
Return all parent communicator topology ranks as np.ndarray,
undistributed axes included.
"""
return np.transpose(
self._topology._proc_ranks_mapping, axes=self._get_proc_axes()
)
def _get_proc_neighbour_ranks(self):
"""
Return the ranks of the neighbours nodes as obtained by MPI_Cart_shift.
self.neighbours[0,i] (resp. [1,i]) is the previous (resp. next) neighbour in axe i.
If axe is not distributed, self.neighbours[:,i] returns [-1,-1].
"""
if self._topology_state.memory_order is MemoryOrdering.F_CONTIGUOUS:
prev, next = self._topology._proc_neighbour_ranks[
:, tuple(self._get_proc_axes())
]
return npw.asintegerarray((prev[::-1], next[::-1]))
else:
return self._topology._proc_neighbour_ranks[:, tuple(self._get_proc_axes())]
def _get_is_distributed(self):
"""
Returns Directions which have been distributed,
is_distributed[dir] = True means that data has been distributed along dir.
"""
return self._proc_transposed(self._topology._is_distributed)
def _get_is_periodic(self):
r"""
MPI communicator grid periodicity.
is_periodic[dir] = True means that the MPI grid is periodic along dir.
/!\ This is not equivalent to domain periodicity, as a periodic
direction might not be distributed in the MPI cartesian grid
or might be forced to be periodic for other reasons through
the is_periodic parameter override.
"""
return self._proc_transposed(self._topology._is_periodic)
def _get_distributed_axes(self):
"""Return distributed axes ids as a np.ndarray of integers."""
return np.where(self._get_is_distributed() == True)[0].astype(np.int32)
def _get_periodic_axes(self):
"""Return cartesian communicator periodic axes ids as a np.ndarray of integers."""
return np.where(self._get_is_periodic() == True)[0].astype(np.int32)
global_resolution = property(_get_global_resolution)
grid_resolution = property(_get_grid_resolution)
ghosts = property(_get_ghosts)
comm = property(_get_comm)
cart_comm = property(_get_cart_comm)
cart_dim = property(_get_cart_dim)
cart_size = property(_get_cart_size)
cart_rank = property(_get_cart_rank)
cart_coords = property(_get_cart_coords)
cart_shape = property(_get_cart_shape)
cart_periods = property(_get_cart_periods)
cart_ranks = property(_get_cart_ranks)
cart_ranks_mapping = property(_get_cart_ranks_mapping)
cart_neighbour_ranks = property(_get_cart_neighbour_ranks)
proc_coords = property(_get_proc_coords)
proc_shape = property(_get_proc_shape)
proc_ranks = property(_get_proc_ranks)
proc_ranks_mapping = property(_get_proc_ranks_mapping)
proc_neighbour_ranks = property(_get_proc_neighbour_ranks)
is_distributed = property(_get_is_distributed)
is_periodic = property(_get_is_periodic)
distributed_axes = property(_get_distributed_axes)
periodic_axes = property(_get_periodic_axes)
[docs]
def default_state(self):
"""Return the default topology state of this topology."""
return CartesianTopologyState(
dim=self.domain.dim, task_id=self.mpi_params.task_id
)
[docs]
def short_description(self):
"""
Returns a short description of the current TopologyView.
Short version of long_description().
"""
s = "{}[domain={}, backend={}, pshape={}, "
s += "grid_resolution={}, ghosts={}, bc=({})]"
s = s.format(
self.full_tag,
self.domain.domain.full_tag,
self.backend.kind,
self.proc_shape,
"[{}]".format(",".join(str(s) for s in self.grid_resolution)),
"[{}]".format(",".join(str(g) for g in self.ghosts)),
",".join(
"{}/{}".format(
str(lb).replace("HOMOGENEOUS_", "")[:3],
str(rb).replace("HOMOGENEOUS_", "")[:3],
)
for (lb, rb) in zip(*self.mesh.global_boundaries)
),
)
return s
[docs]
def long_description(self):
"""
Returns a long description of the current TopologyView.
Long version of short_description().
"""
s = f"======== {self.full_tag} ========\n"
s += " *on task: " + str(self.task_id) + "\n"
s += " *backend: " + str(self.backend.full_tag) + "\n"
s += " *shape: " + str(self.proc_shape) + "\n"
s += " *process of coords " + str(self.proc_coords[:])
s += " and of ranks cart_rank={}, parent_rank={}\n".format(
self.cart_rank, self.mpi_params.rank
)
s += " *cartesian ranks map:\n"
s += prepend(str(self.cart_ranks), " " * 4) + "\n"
s += " *cartesian to parent comm ranks mapping:\n"
s += prepend(str(self.proc_ranks_mapping), " " * 4) + "\n"
s += " *neighbours ranks (left, right) x direction \n"
s += prepend(str(self.proc_neighbour_ranks), " " * 4) + "\n"
s += prepend("*" + str(self.domain), " " * 2)
s += prepend("*" + str(self.mesh), " " * 2)
s += "===================================\n"
return s
[docs]
def can_communicate_with(self, target):
"""True if current topo is complient with target.
Complient means :
- all processes in current are in target
- both topologies belong to the same mpi task
"""
if self.topology == target.topology:
return True
msg = "You try to connect topologies belonging to"
msg += " two different mpi tasks. Set taskids properly or use"
msg += " InterBridge."
assert self.task_id == target.task_id, msg
return self.is_consistent_with(target)
[docs]
def is_consistent_with(self, target):
"""
True if target and current object are equal and
have the same parent. Equal means same mesh, same shape and
same domain.
"""
check_instance(target, CartesianTopologyView)
same_parent = self.parent == target.parent
return npw.equal(self.proc_shape, target.proc_shape).all() and same_parent
[docs]
def __str__(self):
"""Same as self.long_description()."""
return self.long_description()
[docs]
class CartesianTopology(CartesianTopologyView, Topology):
"""
CartesianTopology topologies defined on cartesian meshes which communicates
accross processes through a MPI CartesianTopology communicator.
"""
@debug
def __init__(
self,
domain,
discretization,
mpi_params=None,
cart_dim=None,
cart_shape=None,
is_periodic=None,
cutdirs=None,
mesh=None,
cartesian_topology=None,
cl_env=None,
**kwds,
):
super().__init__(
mpi_params=mpi_params,
domain=domain,
discretization=discretization,
cart_dim=cart_dim,
cart_size=None,
proc_shape=None,
is_periodic=is_periodic,
is_distributed=None,
cartesian_topology=id(cartesian_topology),
mesh=hash(mesh),
topology_state=None,
cl_env=cl_env,
**kwds,
)
[docs]
@debug
def __new__(
cls,
domain,
discretization,
mpi_params=None,
cart_dim=None,
cart_shape=None,
is_periodic=None,
cutdirs=None,
mesh=None,
cartesian_topology=None,
cl_env=None,
**kwds,
):
r"""
Initializes or get an existing CartesianTopology topology.
Parameters
----------
domain : :class:`~hysop.domain.box.Box`
The box geometry on which the cartesian topology is defined.
discretization: :class:`~hysop.tools.parameters.CartesianDiscretization`
Description of the global space discretization of the box (resolution, ghosts,
and boundary conditions).
mpi_params : :class:`~hysop.tools.parameters.MPIParams`, optional
MPI parameters (comm, task ...).
If not specified, comm = domain.task_comm, task = domain.curent_task()
backend: :class:`~hysop.constants.Backend` or `~hysop.core.arrays.ArrayBackend`, optional
Backend or backend kind for this topology.
By default a topology will use Backend.HOST.
cart_dim: int, optional
MPI topology dimension.
cart_shape: list or array of int, optional
MPI grid layout, should be sized as the domain dimension.
is_periodic : tuple, list or array of bool, optional
MPI grid periodicity, *overrides* discretization boundary conditions.
cutdirs: list or array of bool, optional
Set which directions may be distributed,
cutdirs[dir] = True allow MPI to distribute data along dir.
mesh: :class:`~hysop.domain.mesh.CartesianMesh`, optional
A predifined mesh (it includes local and global grids.)
cartesian_topology: :class:`~hysop.core.mpi.MPI.Cartcomm`, optional
A predefined mpi cartesian topology. Use this
when you need the same communicator for two
different meshes/space distribution of data.
kwds: dict
Base class arguments.
Includes allocator, cl_env and queue.
Attributes
----------
global_resolution: np.ndarray of HYSOP_INTEGER
Resolution of the global mesh (as given in the discretization parameter).
ghosts: np.ndarray of HYSOP_INTEGER
CartesianDiscretization ghosts of local-to-process mesh (as given in
the discretization parameter).
mesh: :class:`~hysop.domain.mesh.CartesianMeshView`:
Local mesh on the current mpi process.
cart_comm: :class:`~hysop.core.mpi.MPI.Cartcomm`
MPI cartesian topology intracommunicator.
cart_dim: int
MPI_Cart topology dimension.
Dimension of the MPI cartesian communicator.
This dimension may be smaller than the domain dimension, because
some directions may not be distributed.
cart_size: int
MPI_Cart grid size.
Size of the MPI cartesian communicator (total number of MPI processes).
cart_rank: int
MPI_Cart rank.
Rank of the current process in the cartesian communicator.
May be different than self.mpi_params.rank.
cart_coords: np.ndarray of np.int32
Coordinate of this process in the cartesian communicator grid.
The returned tuple is of dimension self.dimension.
cart_shape: np.ndarray of np.int32
MPI_Cart grid shape.
Shape of the MPI cartesian communicator.
cart_periods: np.ndarray of bool
MPI_Cart grid periodicity
cart_ranks: np.ndarray of np.int32
Return all ranks of this cartesian topology as a np.ndarray such
that array[cart_coords] = rank.
cart_ranks_mapping: np.ndarray of np.int32
Return all ranks of the parent MPI communicator as a np.ndarray such
that array[cart_coords] = parent rank.
cart_neighbour_ranks: np.ndarray
Return the ranks of the neighbours nodes as obtained by MPI_Cart_shift.
self.neighbours[0,i] (resp. [1,i]) is the previous (resp. next)
neighbour in direction i.
proc_coords: tuple of int
Coordinates of this process in the extended cartesian grid
(ie. with non distributed directions included).
The returned tuple is of dimension self.domain_dim.
proc_shape: tuple of int
Processus grid shape, same as cart_shape but extended with non distributed
directions.
proc_ranks: np.ndarray of np.int32
Return all ranks of this cartesian topology as a np.ndarray such
that array[proc_coords] = rank.
proc_ranks_mapping: np.ndarray of np.int32
Return all ranks of the parent MPI communicator as a np.ndarray such
that array[proc_coords] = parent rank.
proc_neighbour_ranks: np.ndarray
Return the ranks of the neighbours nodes as obtained by MPI_Cart_shift.
self.neighbours[0,i] (resp. [1,i]) is the previous (resp. next)
neighbour in axe i.
If axe is not distributed, self.neighbours[:,i] returns [-1,-1].
is_distributed : tuple of bool
Directions which have been distributed,
is_distributed[dir] = True means that data has been distributed along dir.
is_periodic : tuple of bool
MPI grid periodicity.
is_periodic[dir] = True means that the MPI grid is periodic along dir.
/!\ This is not equivalent to domain periodicity, as a periodic
direction might not be distributed in the MPI cartesian grid
or might be forced to be periodic for other reasons through
the is_periodic parameter override.
Domain periodicity is self.domain.periodicity
Notes:
------
* Almost all parameters above are optional.
Only one must be chosen among dim, cutdirs and shape.
See :ref:`topologies` in the Hysop User Guide for details.
* When cartesian_topology is given, dim, shape and cutdirs parameters,
if set, are not used to build the mpi topology, but compared with
cartesian_topology parameters. If they do not fit, error is raised.
* Unless is_periodic is specified periodicity is extracted
from domain boundary conditions.
* All attributes are read-only properties.
"""
# Get or create mpi parameters
mpi_params = cls._create_mpi_params(mpi_params, domain, cl_env)
# prepare MPI processes layout
(cart_dim, cart_size, proc_shape, is_periodic, is_distributed) = (
cls._check_topo_parameters(
mpi_params,
domain,
discretization,
cart_shape,
cutdirs,
cart_dim,
is_periodic,
cartesian_topology,
)
)
# double check types, to be sure RegisteredObject will work as expected
check_instance(mpi_params, MPIParams)
check_instance(domain, Box)
check_instance(discretization, CartesianDiscretization)
check_instance(cart_dim, (int, np.integer))
check_instance(cart_size, (int, np.integer))
check_instance(proc_shape, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(is_periodic, np.ndarray, dtype=bool)
check_instance(is_distributed, np.ndarray, dtype=bool)
check_instance(cartesian_topology, MPI.Cartcomm, allow_none=True)
check_instance(mesh, CartesianMesh, allow_none=True)
cart_dim = int(cart_dim)
cart_size = int(cart_size)
npw.set_readonly(proc_shape, is_periodic, is_distributed)
topology_state = CartesianTopologyState(
dim=domain.dim, task_id=mpi_params.task_id
)
obj = super().__new__(
cls,
mpi_params=mpi_params,
domain=domain,
discretization=discretization,
cart_dim=cart_dim,
cart_size=cart_size,
proc_shape=proc_shape,
is_periodic=is_periodic,
is_distributed=is_distributed,
cartesian_topology=id(cartesian_topology),
mesh=hash(mesh),
topology_state=topology_state,
cl_env=cl_env,
**kwds,
)
if not obj.obj_initialized:
obj.__initialize(
domain,
discretization,
cart_dim,
cart_size,
proc_shape,
is_periodic,
is_distributed,
cartesian_topology,
mesh,
)
return obj
def __initialize(
self,
domain,
discretization,
cart_dim,
cart_size,
proc_shape,
is_periodic,
is_distributed,
cartesian_topology,
mesh,
):
self._discretization = discretization
self._cart_dim = cart_dim
self._cart_size = cart_size
self._proc_shape = proc_shape
self._is_periodic = is_periodic
self._is_distributed = is_distributed
if cartesian_topology is None:
cartesian_topology = self._build_mpi_topo()
self._set_topo(cartesian_topology)
# Get features of the mpi processes grid
self._extract_topo_features(domain)
# Computation of the local meshes
if mesh is None:
mesh = self._compute_mesh(domain, discretization)
self._mesh = mesh
self._TopologyView__set_mesh(mesh)
npw.set_readonly(
self._proc_coords,
self._proc_shape,
self._proc_ranks,
self._proc_neighbour_ranks,
self._is_distributed,
self._is_periodic,
)
if __debug__:
self.__check_vars()
def __check_vars(self):
# check variables and properties at the same time
check_instance(self.global_resolution, np.ndarray, HYSOP_INTEGER)
check_instance(self.ghosts, np.ndarray, HYSOP_INTEGER)
check_instance(self.cart_comm, MPI.Cartcomm)
check_instance(self.cart_dim, int, minval=1)
check_instance(self.cart_size, int, minval=1, maxval=self.mpi_params.size)
check_instance(self.cart_rank, int, minval=0, maxval=self._cart_size)
check_instance(self.cart_coords, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.cart_shape, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.cart_ranks, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.cart_neighbour_ranks, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.cart_ranks_mapping, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.cart_periods, np.ndarray, dtype=bool)
check_instance(self.proc_coords, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.proc_shape, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.proc_ranks, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.proc_neighbour_ranks, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.proc_ranks_mapping, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.distributed_axes, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.periodic_axes, np.ndarray, dtype=HYSOP_INTEGER)
check_instance(self.is_distributed, np.ndarray, dtype=bool)
check_instance(self.is_periodic, np.ndarray, dtype=bool)
check_instance(self.domain, BoxView)
check_instance(self.mesh, CartesianMeshView)
[docs]
def topology_like(
self,
backend=None,
grid_resolution=None,
ghosts=None,
lboundaries=None,
rboundaries=None,
mpi_params=None,
cart_shape=None,
**kwds,
):
"""Return a topology like this object, possibly altered."""
assert "global_resolution" not in kwds, "Specify grid_resolution instead."
grid_resolution = first_not_None(
grid_resolution, self._discretization.grid_resolution
)
ghosts = first_not_None(ghosts, self._discretization.ghosts)
lboundaries = first_not_None(lboundaries, self._discretization.lboundaries)
rboundaries = first_not_None(rboundaries, self._discretization.rboundaries)
backend = first_not_None(backend, self._backend)
discretization = CartesianDiscretization(
resolution=grid_resolution,
ghosts=ghosts,
lboundaries=lboundaries,
rboundaries=rboundaries,
)
# find out the target mpi_params
if __HAS_OPENCL_BACKEND__:
from hysop.core.arrays.all import OpenClArrayBackend
if isinstance(backend, OpenClArrayBackend):
if (mpi_params is not None) and (
mpi_params != backend.cl_env.mpi_params
):
msg = "Backend mpi params mismatch."
raise RuntimeError(msg)
mpi_params = backend.cl_env.mpi_params
mpi_params = first_not_None(mpi_params, self.mpi_params)
if self.domain.dim == grid_resolution.size:
cart_shape = first_not_None(cart_shape, self.proc_shape)
return CartesianTopology(
domain=self._domain,
mpi_params=mpi_params,
discretization=discretization,
backend=backend,
cart_shape=self.proc_shape,
cartesian_topology=None,
**kwds,
)
@classmethod
def _check_topo_parameters(
cls,
mpi_params,
domain,
discretization,
shape,
cutdirs,
dim,
is_periodic,
cartesian_topology,
):
"""
Check compatibility between input parameters provided
to build the mpi topology and return shape of the cartesian topology.
The layout can be obtained by:
(a) - from an already built cartesian topology
(b) - from 'shape' parameter : we choose explicitely the layout.
(c) - from 'cutdirs' parameter: we choose the directions
to split and let MPI fix the number of processes in each dir.
(d) - from dimension of the topology ==> let MPI
choose the 'best' layout.
(e) - in last resort the dimension will be the dimension of the domain.
"""
domain_dim = domain.dim
parent_size = mpi_params.comm.Get_size()
if cartesian_topology:
msg = "Wrong type for input communicator."
assert isinstance(cartesian_topology, MPI.Cartcomm), msg
msg = "The given cartesian topology does not fit with the others"
msg += " input parameters."
assert (dim is None) or (cartesian_topology.dim == dim), msg
assert (shape is None) or (cartesian_topology.size == prod(shape)), msg
assert (shape is None) or (
shape == npw.asintegerarray(cartesian_topology.dims)
).all(), msg
assert is_periodic is not None
dim = cartesian_topology.dim
cart_shape = cartesian_topology.dims
shape = npw.dim_ones(domain_dim)
shape[:dim] = cart_shape
elif shape is not None:
# method (b)
msg = " parameter is useless when shape is provided."
assert cutdirs is None, "cutdirs " + msg
assert dim is None, "dim " + msg
dim = domain_dim
shape = npw.asintegerarray(shape)
msg = "Input shape must be of "
msg += "the same size as the domain dimension."
assert shape.size == dim, msg
elif cutdirs is not None:
# method (c)
msg = " parameter is useless when cutdirs is provided."
assert shape is None, "shape " + msg
assert dim is None, "dim " + msg
shape = npw.dim_ones(domain_dim)
is_distributed = npw.asboolarray(cutdirs).copy()
if is_distributed.any():
dim = np.sum(is_distributed > 0)
assert dim <= domain_dim, "cutdirs is not of size domain_dim"
cart_shape = npw.asintegerarray(MPI.Compute_dims(parent_size, dim))
cls._optimize_shape(cart_shape)
assert (
sum(cutdirs) == cart_shape.size
), "Created shape {} doesnt respect specified cutdirs {}".format(
np.sum(cutdirs > 0), cart_shape.size
)
shape[is_distributed > 0] = cart_shape
else:
assert parent_size == 1
else:
if dim is not None:
# method (d)
msg = " parameter is useless when dim is provided."
assert shape is None, "shape " + msg
assert cutdirs is None, "cutdirs " + msg
else:
# method (e)
dim = domain_dim
cart_shape = npw.asintegerarray(MPI.Compute_dims(parent_size, dim))
shape = npw.dim_ones(domain_dim)
shape[:dim] = cart_shape
cls._optimize_shape(shape)
if is_periodic is None:
try:
is_periodic = discretization.periodicity
except AttributeError:
msg = "Given CartesianDiscretization was not setup correctly:"
msg += "\n => Boundary conditions have not been set."
msg += "\n{}"
msg = msg.format(discretization)
raise ValueError(msg)
shape = npw.asintegerarray(shape)
is_periodic = np.asarray(is_periodic) != 0
assert shape.size == domain_dim
assert is_periodic.size == domain_dim
is_distributed = shape > 1
if parent_size == 1:
is_distributed[-1] = True
cart_shape = shape[is_distributed]
cart_dim = cart_shape.size
cart_size = prod(cart_shape)
is_periodic = is_periodic * is_distributed
assert (cart_dim > 0) and (cart_dim <= domain_dim)
assert prod(cart_shape) == prod(shape)
assert (cart_shape <= shape[is_distributed]).all()
assert cart_size <= parent_size
assert is_periodic.size == domain_dim
assert is_distributed.size == domain_dim
if cart_size < parent_size:
msg = "Only {} processes will be used in the CartesianTopology topology "
msg += "(parent communicator contains {} processes)."
msg = msg.format(cart_size, parent_size)
warnings.warn(msg, TopologyWarning)
proc_shape = shape
return (cart_dim, cart_size, proc_shape, is_periodic, is_distributed)
def _build_mpi_topo(self):
cart_shape = self._proc_shape[self._is_distributed]
periods = self._is_periodic[self._is_distributed]
return self.mpi_params.comm.Create_cart(
cart_shape, periods=periods, reorder=True
)
def _set_topo(self, cartesian_topology):
"""Check and set self._cart_comm and self._cart_rank."""
msg = "Wrong type for input communicator."
assert isinstance(cartesian_topology, MPI.Cartcomm), msg
comm = cartesian_topology
assert self.cart_dim == comm.ndim
assert self.cart_size == comm.Get_size()
assert (
self._proc_shape[self._is_distributed] == npw.asintegerarray(comm.dims)
).all()
assert (
self._is_periodic[self._is_distributed] == npw.asboolarray(comm.periods)
).all()
self._cart_comm = comm
self._cart_rank = comm.Get_rank()
def _extract_topo_features(self, domain):
"""Set self._proc_coords, self._proc_ranks, self._proc_ranks_mapping."""
is_distributed = self._is_distributed
proc_shape = self._proc_shape
cart_comm = self._cart_comm
mpi_params = self._mpi_params
domain_dim = domain.dim
cart_shape = proc_shape[is_distributed]
cart_coords = npw.asintegerarray(cart_comm.coords)
proc_coords = npw.dim_zeros(domain_dim)
proc_coords[is_distributed] = cart_coords
proc_ranks = npw.dim_zeros(proc_shape)
cart_ranks = proc_ranks.reshape(cart_shape)
proc_ranks_mapping = npw.dim_zeros(proc_shape)
cart_ranks_mapping = proc_ranks_mapping.reshape(proc_ranks_mapping.size)
proc_neighbour_ranks = npw.dim_zeros(shape=(2, domain_dim))
for coords in np.ndindex(*cart_shape):
rank = cart_comm.Get_cart_rank(coords)
cart_ranks[coords] = rank
cart_ranks_mapping[...] = MPI.Group.Translate_ranks(
cart_comm.group, cart_ranks.flatten(), mpi_params.comm.group
)
direction = 0
for i in range(self.domain_dim):
if is_distributed[i]:
proc_neighbour_ranks[:, i] = self._cart_comm.Shift(direction, 1)
direction += 1
else:
proc_neighbour_ranks[:, i] = (-1, -1)
assert direction == self.cart_dim
self._proc_coords = npw.asintegerarray(proc_coords)
self._proc_ranks = npw.asintegerarray(proc_ranks)
self._proc_ranks_mapping = npw.asintegerarray(proc_ranks_mapping)
self._proc_neighbour_ranks = npw.asintegerarray(proc_neighbour_ranks)
@staticmethod
def _optimize_shape(shape):
"""
Reorder shape according to the data layout
if arrays are in "fortran" order (column major)
the last dir has priority for distribution.
For C-like (row major) arrays, first dir is the
first to be distributed, and so on.
"""
shape.sort()
if HYSOP_ORDER == "C":
shape[:] = shape[::-1]
def _compute_mesh(self, domain, discretization):
assert isinstance(domain, Box)
assert isinstance(discretization, CartesianDiscretization)
assert (
self.domain_dim == discretization.grid_resolution.size
), "The resolution size differs from the domain dimension."
proc_coords = self._proc_coords
proc_shape = self._proc_shape
# Find out dimension and periodic axes of the domain
domain_dim = domain.dim
periodicity = discretization.periodicity
# /!\ Now we assume that the user gives us the grid resolutionn
# and not the global_resolution as it used to be.
# We do not remove 1 point on each periodic axe because of periodicity
computational_grid_resolution = discretization.grid_resolution
# Number of "computed" points (i.e. excluding ghosts).
# /!\ we try to match fftw_mpi_local_size_* functions for the Fortran spectral backend
assert all(computational_grid_resolution >= proc_shape)
pts_noghost = npw.dim_zeros(domain_dim)
pts_noghost[:] = npw.ceil(
npw.divide(computational_grid_resolution.astype(HYSOP_REAL), proc_shape)
)
# If any, remaining points are added on the mesh of the last process.
remaining_points = npw.dim_zeros(domain_dim)
remaining_points[:] = (
computational_grid_resolution - (proc_shape - 1) * pts_noghost
)
assert (remaining_points >= 1).all(), remaining_points
# Total number of points (size of arrays to be allocated)
nbpoints = pts_noghost.copy()
for i in range(domain_dim):
if proc_coords[i] == proc_shape[i] - 1:
nbpoints[i] = remaining_points[i]
local_resolution = nbpoints.copy()
local_resolution += 2 * discretization.ghosts
msg = "\nLocal compute shape is smaller than the total number of ghosts, "
msg += "on at least one axis:"
msg += f"\n *compute shape: {nbpoints}"
msg += f"\n *ghosts: 2*{discretization.ghosts}"
if not np.all(nbpoints >= 2 * discretization.ghosts):
raise RuntimeError(msg)
# Global indices for the local mesh points
global_start = proc_coords * pts_noghost
return CartesianMesh(
topology=self, local_resolution=local_resolution, global_start=global_start
)
[docs]
def view(self, topology_state):
"""
Returns a view of this topology with the given state.
"""
check_instance(topology_state, CartesianTopologyState)
return CartesianTopologyView(topology=self, topology_state=topology_state)
[docs]
def discretize(self, field):
"""Discretize a continous field on this topology and return a DiscreteField."""
from hysop.fields.continuous_field import ScalarField
from hysop.fields.cartesian_discrete_field import (
CartesianDiscreteScalarField,
TmpCartesianDiscreteScalarField,
)
check_instance(field, ScalarField)
if (field.lboundaries_kind != self._discretization.lboundaries).any() or (
field.rboundaries_kind != self._discretization.rboundaries
).any():
msg = """
Cannot discretize a field with cartesian boundary conditions:'
lboundaries: {}
rboundaries: {}
On a cartesian topology with different boundary conditions:
lboundaries: {}
rboundaries: {}
""".format(
field.lboundaries_kind,
field.rboundaries_kind,
self._discretization.lboundaries,
self._discretization.rboundaries,
)
raise RuntimeError(msg)
if field.is_tmp:
return TmpCartesianDiscreteScalarField(field=field, topology=self)
else:
return CartesianDiscreteScalarField(field=field, topology=self)
